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We perform 3D numerical simulations in full general relativity to study the stability of rapidly 
rotating, supramassive neutron stars at the mass-shedding limit to dynamical collapse. We adopt 
an adiabatic equation of state with T = 2 and focus on uniformly rotating stars. We find that 
the onset of dynamical instability along mass-shedding sequences nearly coincides with the onset of 
secular instability. Unstable stars collapse to rotating black holes within about one rotation period. 
We also study the collapse of stable stars which have been destabilized by pressure depletion (e.g. via 
a phase transition) or mass accretion. In no case do we find evidence for the formation of massive 
disks or any ejecta around the newly formed Kerr black holes, even though the progenitors are 
rapidly rotating. 
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I. INTRODUCTION 

Neutron stars found in nature are rotating. Rotation 
can support stars with higher mass than the maximum 
static limit, producing "supramassive" stars, as defined 
and calculated by Cook, Shapiro and Teukolsky 
Such supramassive stars can be created when neutron 
stars accrete gas from a normal binary companion. This 
scenario can also lead to "recycled" pulsars [see || for 
model calculations in general relativity]. Alternatively, 
supramassive stars can be produced in the merger of bi- 
nary neutron stars ( || for discussions and references). 

Pulsars are believed to be uniformly rotating. Even- 
tually, viscosity will drive any equilibrium star to uni- 
form rotation. Uniformly rotating configurations with 
sufficient angular momentum will be driven to the mass- 
shedding limit (at which the star's equator rotates with 
the Kepler frequency, so that any further spin-up would 
disrupt the star; ||). Supramassive neutron stars at the 
mass-shedding limit is the subject of this paper. 

The dynamical stability of rotating neutron stars, in- 
cluding supramassive configurations, against radial per- 
turbations, as well as the final fate of unstable stars un- 
dergoing collapse, has not been established definitively. 



Along a sequence of nonrotating, spherical stars, pa- 
rameterized by central density, the maximum mass con- 
figuration defines a critical density above which the stars 
are unstable to radial oscillations: stars on the high den- 
sity, unstable branch collapse to black holes on dynamical 
timescales || f?J. 

Establishing the onset of instability for rotating stars is 
more complicated. Chandrasekhar and Friedman @] and 
Schutz |j] developed a formalism to identify points of dy- 
namical instability to axisymmetric perturbations along 
sequences of rotating stars. In their formalism, however, 
a complicated functional for a set of trial functions has 
to be evaluated. Probably because of the complexity of 
their methods, explicit calculations have never been per- 
formed. Friedman, Ipser and Sorkin [Io| showed that 
for uniformly rotating stars the onset of secular instabil- 
ity can be located quite easily by applying turning-point 
methods along sequences of constant angular momentum. 
This method has been applied to find points of secular 
instability in numerical models of uniformly rotating neu- 
tron stars 0^0. 

Turning point methods along such sequences can only 
identify the point of secular, and not dynamical insta- 
bility, since one is comparing neighboring, uniformly ro- 
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tating configurations with the same angular momentum. 
Maintaining uniform rotation during perturbations tac- 
itly assumes high viscosity. In reality, the star will pre- 
serve circulation as well as angular momentum in a dy- 
namical perturbation, and not uniform rotation. It is 
thus possible that a secularly unstable star may be dy- 
namically stable: for sufficiently small viscosity, the star 
may change to a differentially rotating, stable configura- 
tion instead of collapsing to a black hole. Ultimately, the 
presence of viscosity will bring the star back into rigid 
rotation, driving the star to an unstable state. A secu- 
lar instability evolves on a dissipative, viscous timescale, 
while a dynamical instability evolves on a collapse (free- 
fall) timescale. Friedman, Ipser and Sorkin |lCj ] showed 
that along a sequence of uniformly rotating stars, a sec- 
ular instability always occurs before a dynamical insta- 
bility (implying that all secularly stable stars are also 
dynamically stable). 

For spherical stars, the onset of secular and dynamical 
instability coincides (since for a nonrotating star a radial 
perturbation conserves both circulation and uniform ro- 
tation). This suggests that for uniformly rotating stars 
for which the rotational kinetic energy T is typically a 
small fraction of the gravitational energy W, the onset 
of dynamical instability is close to the onset of secular 
instability. One goal of this paper is to test this hypoth- 
esis and to identify the onset of dynamical instability in 
rotating stars. 

We also follow the nonlinear growth of the radial insta- 
bility and determine the final fate of unstable configura- 
tions. Nonrotating neutron stars collapse to black holes, 
but rotating stars could also form black holes surrounded 
by massive disks. Also, if J/Mg exceeds the Kerr limit 
of unity (where J is the angular momentum and M g the 
total mass-energy or gravitational mass of the progeni- 
tor star), not all of the matter can collapse directly to a 
black hole. As pointed out recently, such a system could 
be the central source of 7-ray bursts [H . 



Numerical hydrodynamic simulations in full general 
relativity (GR) provide the best approach to under- 
standing the collapse of rotating neutron stars. Two 
groups |l5| , |l4| ] included rotation in axisymmetric, rela- 
tivistic hydrodynamic codes to study the collapse of ro- 
tating massive stars to black holes. The collapse and 
fate of unstable rotating neutron stars, however, has 
never been simulated before. Probably this is because 
numerical methods for constructing initial data describ- 
ing rapidly rotating neutron stars, as well as numerical 
tools, techniques and sufficient computational resources 
have only become available recently. Over the last few 
years, robust numerical techniques for constructing equi- 
librium models of rotating neutron stars in full GR have 
been established ]i|,p|,pT 15-1^]. More recently, methods 
for the numerical evolution of 3D gravitational fields have 
been developed (see, e.g., ]l8|-p4|]). In a previous pa- 
per p5[ , Shibata presented a wide variety of numerical 
results of test problems for his 3D hydrodynamic GR 
code and demonstrated that simulations for many inter- 
esting problems are now feasible. 

In this paper, we perform simulations in full GR for 
rapidly rotating neutron stars. This study is a by- 
product of our long-term effort to build robust, fully rel- 
ativistic, hydrodynamic evolution codes in 3D. We adopt 
rigidly rotating supramassive neutron stars at mass- 
shedding as initial data. By exploring rotating stars 
at mass-shedding, we can clarify the effect of rotation 
most efficiently. Such stars are also the plausible out- 
come of pulsar recycling and binary coalescence. Fol- 
lowing Ref. p^| , we prepare equilibrium states for such 
stars using an approximation in which the spatial metric 
is assumed to be conformally flat. We then perform nu- 
merical simulations to investigate the dynamical stability 
of the rapidly rotating neutron stars against collapse and 
to determine the final fate of the unstable neutron stars. 
We believe that this is the first 3D simulation of the dy- 
namical collapse of a rotating neutron star in full GR. 
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In Newtonian physics, stars with sufficient rotation 
(T/|W| > 0.27) are dynamically unstable to bar forma- 
tion [p6|-p8|. Since T/|W| increases approximately with 
i?" 1 as a star collapses, radial collapse may drive the 
dynamical growth of nonaxisymmetric bars. To allow 
for this possibility, a numerical simulation must be per- 
formed in 3D, not in axisymmetry. 

In Sec. II, we briefly describe our formulation, ini- 
tial data, and spatial gauge conditions. In Sec. Ill, we 
present numerical results. First, we study the dynamical 
stability of supramassive rotating neutron stars at the 
mass-shedding limit. We then study the final products 
of the unstable neutron stars adopting three kinds of ini- 
tial conditions: In the first case, we choose a marginally 
stable neutron star and slightly reduce the pressure for 
destabilization as the initial condition. In the second 
case, we prepare a stable star, and then reduce a large 
fraction of the pressure suddenly. While we are primarily 
interested in computational consequences, this scenario 
may provide a model for sudden phase transitions inside 
neutron stars (see, e.g., |2^,Q and references therein). 
In the third case, we prepare a stable star and add more 
mass near the surface to induce collapse. In all the cases, 
we find that the final products are black holes without 
surrounding massive disks, which we can readily explain. 
In Sec. IV we provide a summary. Throughout this pa- 
per, we adopt the units G = c = Mq = 1 where G, 
c and Mq denote the gravitational constant, speed of 
light and solar mass, respectively. We use Cartesian co- 
ordinates x k = (x,y,z) as the spatial coordinate with 
r = \J x' 2 + y 2 + z 2 -, t denotes coordinate time. 

II. METHODS 
A. Summary of formulation 

We perform numerical simulations using the same for- 
mulation as in |25|| , to which the reader may refer for 



details and basic equations. The fundamental variables 
used in this paper are: 

p: rest mass density, 

e: specific internal energy, 
P: pressure, 
m m : four velocity, 

a: lapse function, 
[3 k : shift vector, 

jij: metric in 3D spatial hypersurface, 
7= e 120 = det( 7y ), 

lij = e Tij i 

Kij \ extrinsic curvature. 

General relativistic hydrodynamic equations are solved 
using the van Leer scheme for the advection terms fl3l| . 
Geometric variables (together with three auxiliary func- 
tions Fi and the trace of the extrinsic curvature) are 
evolved with a free evolution code. The boundary con- 
ditions for geometric variables are the same as those 
adopted in pl| . Violations of the Hamiltonian constraint 
and conservation of mass and angular momentum are 
monitored as code checks. Several test calculations, in- 
cluding spherical collapse of dust, stability of spherical 
neutron stars, and the evolution of rotating neutron stars 
as well as corotating binary systems have been described 
in |25|] . Black holes that form during the late phase of 
the collapse are located with an apparent horizon finder 
as described in |$2| . 

We also define a density p*(= pau°e 6 ^) from which the 
total rest mass of the system can be integrated as 

M* = J d 3 xp*. (2.1) 

We perform the simulations assuming 7r-rotation sym- 
metry around the z-axis as well as a reflection symmetry 
with respect to the z = plane and using a fixed uni- 
form grid with the typical size 153 x 77 x 77 in x — y — z. 
We have also performed test simulations with different 
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grid resolutions to check that the results do not change 
significantly. 

The slicing and spatial gauge conditions we use in this 
paper are basically the same as those adopted in p5[ ; 
i.e., we impose an "approximate" maximal slice condition 
(K k k ~ 0) and an "approximate" minimum distortion 
gauge condition (D i (d t j l: ') ~ where £>, is the covariant 
derivative with respect to 7y). However, for the case 
when a rotating star collapses to a black hole, we slightly 
modify the spatial gauge condition in order to improve 
the spatial resolution around the black hole. The method 
of the modification is described in Sec. II. C. 

B. Initial conditions for rotating neutron stars 

As initial conditions, we adopt rapidly and rigidly 
rotating supramassive neutron stars in (approximately) 
equilibrium states. The approximate equilibrium states 
are obtained by choosing a conformally flat spatial met- 
ric, i.e., assuming jij = e 4 ^8ij. This approach is compu- 
tationally convenient and, as illustrated in ]35|], provides 
an excellent approximation to exact axisymmetric equi- 
librium configurations. 

Throughout this paper, we assume a T-law equation of 
state in the form 



P = (T - l)pe, 



(2.2) 



where T is the adiabatic constant. For hydrostatic prob- 
lems, the equation of state can be rewritten in the poly- 
tropic form 

P = Kp r , r=l + ~, (2.3) 
n 

where K is the polytropic constant and n the polytropic 
index. We adopt T = 2 (n = I) as a reasonable quali- 
tative approximation to realistic (moderately stiff) cold, 
nuclear equations of state. 

Physical units enter the problem only through the 
polytropic constant K, which can be chosen arbitrarily 



or else completely scaled out of the problem. We of- 
ten quote values for K = 200/7T, for which in our units 
(G=c = M Q = I) the radius is R = (nK/2) 1 / 2 = 10 in 
the Newtonian limit; corresponding to R ~ 15 km. Since 
K n / 2 has units of length, dimensionless variables can be 
constructed as 



M* = M*K-' l/2 ,M g = M g K- n / 2 , R = RK~ n/2 , 

J = JK- n , P = PJf-"/ 2 , p = pK n , (2.4) 



where P denotes rotational period. All results can be 



scaled for arbitrary K using Eqs. (2.4). 

For the construction of the (approximate) equilibrium 
states as initial data, we adopt a grid in which the semi 
major axes of the stars, along the x and y-ax.es, are cov- 
ered with 40 grid points. For rotating stars at mass- 
shedding near the maximum mass, the semi minor (rota- 
tion) axis along the z-axis is covered with 23 or 24 grid 
points. 

In Fig. I, we show the relation between the grav- 
itational mass M g and the central density p c for the 
neutron stars. The solid and dotted lines denote the 
relations for spherical neutron stars and stars rotating 
at the mass-shedding limit, constructed from the exact 
stationary matter and field equations. The open cir- 
cles denote the approximate equilibrium states at the 
mass-shedding limit obtained using the conformal flat- 
ness approximation. We find that at p c — p max where 
0.0040 < Pmax 0.0045, M g takes its maximum value. 
For stars with stiff equations of state the numerical re- 
sults in Ref. show that the central density at the onset 
of secular instability (which hereafter we refer to as p cr it) 
is very close to p max (see, e.g., Fig. 4 of j|). We therefore 
consider stars with p c > p max (~ p C rit) as candidates for 
dynamically unstable stars. 
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C. Spatial gauge condition 

When no black hole is formed, we adopt the approx- 
imate minimum distortion gauge condition as our spa- 
tial gauge condition (henceforth referred to as the AMD 
gauge condition). However, as pointed out in previous 
papers p^ , p5| , during the black hole formation (i.e., for 
an infalling radial velocity field), the expansion of the 
shift vector dij3 l and the time derivative of <fr becomes 
positive using this gauge condition together with max- 
imal slicing. Accordingly, the coordinates diverge out- 
ward and the resolution around the black hole form- 
ing region becomes worse and worse during the collapse. 
This undesirable property motivates us to modify the 
AMD gauge condition when we treat black hole forma- 
tion. Specifically, we modify the AMD shift vector ac- 
cording to 



ft 1 — /^amd - /(*> r ) — ; — /^amd- 
r + e 



(2.5) 



Here /3X MD denotes the shift vector in the AMD gauge 
condition, /?Xmd = xk Pamd / ( r ~^ £ ) ' e is a small constant 
much less than the grid spacing, and f(t, r) is a function 
chosen as 



f(t,r) = / (i): 



1 



(2.6) 



l + (r/M ff , ) 4 ' 
where M g $ denotes the gravitational mass of a system 
at t = 0. We determine fo(t) from (fro — <fr(r = 0). 
Taking into account the fact that the resolution around 
r = deteriorates when (fro becomes large, we choose fo 
according to 

1 for (fr > 0.8, 

fo{t) = { 2.50o - 1 for 0.4 < cfr < 0.8, Type I (2.7) 
for (fro < 0.4, 



1 for (fro > 0.6, 

fo(t) = { 5</>o-2 for 0.4 < <fr < 0.6, Type II (2.8) 
for (fro < 0.4. 

Note that for spherical collapse with /q = 1, di/3 % ~ at 
r = in both cases. In general, we find numerically that 



dif} % is small near the origin, where the collapse proceeds 
nearly spherically. In the following, we refer to the modi- 
fied gauge conditions of /o defined by Eqs. ( |2.7|) and (2.S) 



as type I and II, respectively. We employ them whenever 
a rotating neutron star collapses to a black hole. 

III. NUMERICAL RESULTS 
A. Dynamical stability 

We investigate the stability of supramassive rotating 
neutron stars at mass-shedding limits against gravita- 
tional collapse. We adopt the stars marked with (A), 
(B), (C), (D), and (E) in Fig. 1 as initial conditions for 
our numerical experiments. The physical properties of 
these stars are summarized in Table I. In this numeri- 
cal experiment, we adopt two initial conditions for each 
model. In one case, we use the (approximate) equilibrium 
states of rotating neutron stars without any perturbation 
and in the other case, we uniformly reduce the pressure 
by 1% (by decreasing K; i.e., AK/K = 1% where AK 



denotes the depletion factor of K 35 1). 

In Fig. 2, we show p and a at r = |34| as a function of 
t/P where P is the rotation period of each rotating star. 
We find that when p c < p crit (i.e., stars (A), (B) and 
(C)), the rotating stars oscillate independent of the ini- 
tial perturbations. Hence, these stars are stable against 
gravitational collapse. We note that we find small ampli- 
tude oscillations even when we do not reduce the pressure 
initially. This is caused by small deviations of the initial 
data from true equilibrium states, both because of the 
conformal flatness approximation and because of numer- 
ical truncation error. 

We expect the oscillation frequencies in Fig. 2 to be 
the fundamental quasi-radial oscillation of these rotating 
stars. The oscillation periods increase with the central 
density. At the marginally stable point of secular stabil- 
ity (p c = Pcrit), the period becomes infinite. 
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Star (D) does not collapse either and instead oscillates 
for AK = 0. However, it collapses for AK/K = 1%. 
This indicates that star (D) is located near the onset 
point of dynamical stability. It is found that the oscilla- 
tion amplitude for the case AK = is very large com- 
pared with those for (A)-(C). This could be caused by 
two effects: (i) star (D) is near the onset of dynamical 
instability and hence a small deviation from true equi- 
librium can induce large perturbations; (ii) the confor- 
mal flatness approximation results in larger deviations 
from true equilibrium for more relativistic configurations, 
which causes a larger initial perturbation. Apparently, 
the deviation caused by the numerical truncation error 
and/or the conformal flatness approximation stabilizes 
the configuration, and for AK = the star oscillates with 
an average value of the central density (p(r = 0) ~ 0.004) 
slightly smaller than the initial value p c ~ 0.0047. This 
suggests that star (D) with AK = is a perturbed state 
of a true equilibrium star of p c ~ 0.004 ~ p cr it. The 
results for star (E) are similar to, but more pronounced 
than those for star (D), suggesting that the initial con- 
figuration (E) may also be a perturbation of a stable star 
with p c ~ 0.004 ~ p cr it. 

To determine the onset of dynamical instability 
more sharply, we perform further simulations adopting 
AK/K = 0.7%, 0.8%, and 0.9% for star (D). In Fig. 3, 
we show p and a at r = as a function of t/P for star (D) 
for the various initial depletion factors. We find that for 
AK/K < 0.8%, the stars behave similarly to AK = 0; 
i.e., the stars simply oscillate with the average density 
— Pcrit- For AK/K > 0.9%, however, the stars quickly 
collapse to a black hole. We do not find any examples 
in which the stars oscillate with average densities larger 
than 0.0045 > Pcrit- This indicates that equilibrium stars 
with p c > /9 cr ;t are dynamically unstable. Although we 
cannot specify the onset of dynamical instability with 
arbitrary precision, our present results indicate that it 
nearly coincides with the onset of secular instability. 



B. Final outcome of unstable collapse 

To study the final outcome of the gravitational collapse 
of rapidly rotating neutron stars, we perform a number 
of numerical experiments for several different initial con- 
ditions. 

First, we evolve star (D) with AK/K = 1% for differ- 
ent spatial gauge conditions. In Fig. 4, we show </> and a 
at r = as a function of time for the AMD gauge (solid 
line), the modified gauge of type I (dotted line) and type 
II (dashed line). As stated in Sec. II. B, <p(r = 0) in- 
creases quickly during the gravitational collapse for the 
AMD gauge. In this case, a(r = 0) stops decreasing in 
the late phase of the collapse where <p(r = 0) > 1, which 
is a numerical artifact. This is probably caused by the in- 
sufficient resolution around the black hole forming region. 
For the modified gauge conditions, a(r = 0) smoothly 
approaches zero. We note that a(r = 0) ought to be 
independent of the spatial gauge condition, so that the 
deviation of the AMD results from the modified gauge 
condition results are a numerical artifact. This shows 
that the results for t/P > 1.4 computed in the AMD 
gauge condition is unreliable and indicates that the mod- 
ification of the AMD gauge condition is an appropriate 
strategy to overcome the deterioration of the resolution 
in the late phase of the collapse. 

In Fig. 5, we show the time variation of the total angu- 
lar momentum of the system. Since the evolving system 
is nearly axisymmetric, the angular momentum should 
be nearly conserved. In all the three cases, however, the 
angular momentum slowly decreases in the early phase, 
which is caused by numerical dissipation at the stellar 
surface. As the collapsing star approaches a black hole, 
the angular momentum changes quickly because the res- 
olution becomes increasingly worse. In the AMD gauge 
case, the error amounts to > 5%, while in the modified 
gauge cases, it is ~ 1.5% at the time when apparent 
horizon is found at t ~ 1.4P (see Fig. 6). This is further 
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evidence that the modified gauge conditions are better 
suited for simulations of black hole formation. 

It should be noted that even with the modified gauge 
conditions, the resolution becomes too poor to perform 
accurate simulations for times exceeding t/P ^ 1.5. This 
is because the metric 7^ becomes very spiky around the 
apparent horizon (i.e., because of horizon throat stretch- 
ing). To perform simulations for times much later than 
horizon formation special computational tools are neces- 
sary, probably including apparent horizon boundary con- 
ditions n). 

In Fig. 6, we show snapshots of density contour lines 
for the density /?* and the velocity field for v l (= u l /u°) in 
the equatorial and y — planes. The results are obtained 
in the modified gauge condition of type I. It is found that 
after about 1.4 orbital periods almost all the matter has 
collapsed to a black hole. In Fig. 7, we show the fraction 
of the rest mass inside a coordinate radius r, defined as 

M*(r) 



1 



d xp* 



(3.1) 



'* J\x'\<r 

R e denotes the coordinate axial length in the equatorial 
plane at t = (see Table I). Note that at t ~ 1.4P, 
the apparent horizon is located at r ~ 0.2R e . Thus, 
almost all the matter (more than 99%) has been ab- 
sorbed by the black hole by that time. Although Fig. 
6 shows that a small fraction of the matter has not yet 
been swallowed by the black hole, the matter which stays 
inside r < R e ~ 5M g will ultimately have to fall in. 
This is, because the radius of the innermost stable cir- 
cular orbit (ISCO) is i?fsco ~ ^M g for a (nonrotating) 
Schwarzschild black hole in our gauge. The collapse of 
rotating neutron stars with J/Mg ~ 0.6 < 1 leads to 
moderately rotating Kerr black holes, for which Rfg CO 
is an adequate approximation to the ISCO. This fact al- 
ready suggests that no disk will form around the black 
hole. 

The same reason also suggests why no massive disk 
forms during the collapse: the equatorial radius R e is 



initially less than 5M ff , and hence inside the radius which 
will become the ISCO of the final black hole. 

Next, we evolve the initial configuration (A) deplet- 
ing the pressure by various amounts, which may pro- 
vide a model for sudden phase transitions inside neutron 
stars |2l||3"(i| ] . In Fig. 8, we show p and a at r = as a 
function of time for AK/K = 0, 1%, 5%, and 10% @. 
When the depletion factor is less than 5%, the star sim- 
ply oscillates, but for AK/K = 10% the star collapses 
dynamically. Note that depleting the pressure by 10% is 
approximately equivalent to increasing the gravitational 
mass by 5% according to the scaling relation for the poly- 



tropic stars of r = 2 (see Eq. (2.4)). Since the gravita- 
tional mass for star (A) is about 3% less than the max- 
imum allowed mass, it is quite reasonable that this star 
collapses. In the following two simulations, we focus on 
evolutions of star (A) with AK/K = 10%. 

In order to test if nonaxisymmetric (bar-mode) pertur- 
bations have enough time to grow appreciably during the 
gravitational collapse, we excite such a perturbation by 



modifying the initial density profile p* according to 35 

. x — y 



P* = (p*)o 1 + 0.3 



(3.2) 



where (p*)o denotes the density profile of star (A) in the 
unperturbed state. 

In Figs. 9 and 10, we show snapshots of density contour 
lines for and the velocity field for v % in the equato- 
rial and y = planes for the above axisymmetric and 
nonaxisymmetric initial conditions. For these simula- 
tions we adopted the modified gauge condition of type 
II. In Fig. 11, we also show M*(r)/M* as a function of 
time for these cases. We again find that irrespective of 
the initial perturbation, almost all the matter collapses 
into the black hole without any massive disk or ejecta 
around the black hole. Again, this is a consequence of 
the stars being sufficiently compact that almost all the 
matter ends up inside the ISCO of the final black hole. 
Note that the star with the nonaxisymmetric perturba- 
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tion evolves very similarly to the unperturbed, axisym- 
metric star, showing that the dynamical collapse does not 
leave the nonaxisymmetric perturbation enough time to 
grow appreciably during the collapse. Again, this can be 
understood quite easily from the following heuristic (and 
Newtonian) argument. Star (A) has an initial equatorial 
radius of R e ~ 5.5M S , and can therefore shrink by less 
than a factor of 3 before a black hole forms. Its initial 
value of T/\W\ is about 0.09 (see Table I). Since T/\W\ 
scales approximately with R , it can just barely reach 
the critical value (T/|W|)dyn ~ 0.27 for dynamical in- 
stability before a black hole forms. It is therefore not 
surprising that we do not find dynamically growing ax- 
isymmetric perturbations. Note that the star does reach 
the critical value for secular instability to bar formation 
(which may be as small as (T/|W|) sec ~ 0.1 for very com- 
pact configurations, see |37j), so that viscosity or emis- 
sion of gravitational waves could drive the star unstable. 
However, this mode would grow on the corresponding 
dissipative timescale, which is much longer than the dy- 
namical timescale of the collapse. 

In order to make these statements about nonaxisym- 
metric growth more quantitative, we compare the quan- 
tities 



p* 



(p*)o( 



1 + 0.5 



(3.5) 



^rms ~r l/rms 



(3-3) 



i ^2 



drxp*(x l ) 



1/2 



(3.4) 



for the perturbed and unperturbed evolutions in Fig. 12. 
Here, x\ ms denotes the mean square axial length defined 

as 

1 

The figure shows very clearly that the axial ratio oscil- 
lates for the perturbed evolution, but does not grow on 
the dynamical timescale of the collapse. 

Finally, we model a scenario in which a small amount 
of matter accretes onto a stable star resulting in desta- 
bilization of the star. As the stable star, we again adopt 
configuration (A) and to model the matter accretion we 
modify the initial density distribution according to |35| 



with all the matter moving with the same initial angular 
velocity. Most of the enhancement is in the outer region, 
which mimics the effect of accretion. In this case, the 
total rest mass is about 9.5% larger than that of star 
(A), so that the mass is larger than the maximum allowed 
mass along the sequence of rotating neutron stars. The 
value of J/Mg is nearly unchanged. Note that we do 
not reduce the pressure for these simulations. We again 
evolve the star using the modified gauge condition of type 
II. 

The star again evolves very similarly to those in the 
previous two cases. As an example, we show in Fig 13 
M*(r)/-M* as a function of time, which is similar to the 
results in Fig. 11. The apparent horizon forms at t ~ 
0.89P and r ~ 0.2R e . We again find that almost all the 
matter collapses into the black hole without forming a 
massive disk around the black hole. 

IV. SUMMARY AND CONCLUSION 

We perform fully relativistic, 3D hydrodynamic simu- 
lations of supramassive neutron stars rigidly rotating at 
the mass-shedding limit. We study the dynamical stabil- 
ity of such stars close to the onset of secular instability 
and follow the collapse to rotating black holes. 

Our results suggest that the onset of dynamical, radial 
instability is indeed close to the onset of secular instabil- 
ity, as expected from the coincidence of the secular and 
dynamical instability in nonrotating, spherical stars. 

In all our simulations, nearly all the matter is con- 
sumed by the nascent black hole by the time the calcula- 
tion stops, and we do not find any evidence for a forma- 
tion of a massive disk or any ejecta. Since we are con- 
sidering maximally rotating neutron stars at the mass- 
shedding limit, and since the formation of a disk is even 
less likely for more slowly rotating stars, we conclude that 
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such disks quite generally do not form during the collapse 
of unstable, uniformly rotating neutron stars. This also 
includes stars which are destabilized by pressure deple- 
tion (as, for example, by a nuclear phase transition), or 
by mass accretion. 

We also find that during the collapse to a black hole, 
nonaxisymmetric perturbations do not have enough time 
to grow appreciably. 

Both these findings can be understood quite easily 
from heuristic arguments. The initial equilibrium con- 
figurations are sufficiently compact, typically R e < 6M g , 
so that most of the matter already starts out inside the 
radius which will become the ISCO of the final black 
hole. Therefore it is very unlikely that a stable, mas- 
sive disk would form. Also, the star can only contract 
by about a factor of three before a black hole forms. 
Hence T/|W|, which approximately scales with can 
only increase by about a factor of three over its initial 
value of (Ty|W|)init ~ 0.09, and only barely reaches the 
critical value of dynamical instability for bar formation 
(r/|W|)dy n ^0.27. It is therefore not surprising that 
we do not see a dynamical growth of nonaxisymmetric 
perturbations. We expect that these results hold for any 
moderately stiff equation of state, for which the corre- 
sponding critical configurations are similarly compact. 

The study reported here focuses on uniformly rotat- 
ing neutron stars, for which we adopt a moderately stiff 
equation of state and consider a configuration which is 
moderately compact initially (R/M g ~ 6). We specu- 
late that for two alternative scenarios the results may 
be quite different, even qualitatively, both as far as the 
formation of a disk and the growth of nonaxisymmetric 
perturbations are concerned. 

For differentially rotating neutron stars, which are the 
likely outcome of the merger of binary neutron stars [Q , 
T/|W| may take larger values than for rigidly rotating 
neutron stars. It is therefore possible that such stars 
might develop dynamical bar mode instabilities. 



Rotating supermassive stars (with masses M > 
10 5 Mq) or massive stars on the verge of supernova col- 
lapse are subject to the same dynamical instabilities, 
but are characterized by very soft equations of state 
(r ~ 4/3) and initial configurations which are nearly 
Newtonian (see |7j,[38|]). Such stars therefore reach the 
critical value (T/|W|)d yn for bar mode formation far out- 
side the horizon radius. Moreover, R/M is very large 
initially, so that a disk may easily form (compare the 
discussion in @). 

We will treat the collapse of both differentially rotating 
neutron stars and supermassive stars in future papers. 
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Table I. The list of the central density p c , total rest 
mass M*, gravitational mass M g , angular momentum in 
units of (J/Mg), T/W, rotation period P, and coor- 
dinate length of the semi-major axis R e for rotating neu- 
tron stars at mass-shedding limits in the conformal flat 
approximation. The gauge invariant definition of T/W 
is the same as that in Ref. jlj. The units of mass, length 
and time are M Q , 1.477km, and 4.927^isec, respectively. 



Pc(l(T J ) 


M* 


M g 




T/W 


P 


Re 


Model 


2.77 


1.580 


1.452 


0.598 


0.087 


163 


8.064 


(A) 


3.38 


1.628 


1.484 


0.586 


0.085 


148 


7.820 


(B) 


3.98 


1.646 


1.496 


0.574 


0.083 


137 


7.365 


(C) 


4.68 


1.645 


1.494 


0.563 


0.080 


127 


6.934 


(D) 


5.43 


1.630 


1.481 


0.553 


0.078 


120 


6.566 


(E) 




0.006 



FIG. 1. The gravitational mass M g as a function of 
the central density p c for rotating stars with T = 2 and 
K = 200/-7T. The solid and dashed lines denote exact solutions 
for sequences of rotating stars at the mass-shedding limit and 
spherical stars. The open circles denote the sequence of ro- 
tating stars at the mass-shedding limit as obtained from the 
conformal flatness approximation. The configurations that we 
adopt as initial data for dynamical evolution calculations in 
this paper are marked with (A)-(E). 
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FIG. 2. a and p at r = as a function of t/P in the evolu- 
tion of stars (A)-(E). The solid and dotted lines denote the 
results for AK = and AK/K — 1%, respectively. 



FIG. 4. <f> and a at r — as a function of t/P during the 
collapse of star (D) with AK/K = 1% for the AMD gauge 
(the solid line) , the modified gauge of type I (the dotted line) 
and of type II (the dashed line), respectively. 
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FIG. 3. a and p at r = as a function of t/P during the 
evolution of star (D) for various AK/K. The solid, dashed, 
long dashed, dotted-dashed, and dotted lines denote the re- 
sults for AK/K = 0, 0.7%, 0.8%, 0.9% and 1%. 



FIG. 5. Same as Fig. 4, but for the angular momentum 
J / Jo as a function of t/P. Here, Jo is the angular momentum 
of the system at t — 0. 
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t = O.OOP 




FIG. 6. Snapshots of density contours for p» and the ve- 
locity flow for (v x ,v v ) in the equatorial plane (left) and 
in the y = plane (right) for the collapse of star (D) 
with AK/K = 1% (evolved with type I modified AMD 
gauge). The contour lines are drawn for p*/p» c = lO -0 ' 3 ^ 
for j = 0, 1, 2, • • • , 10 where p, c is 0.034, 0.64 and 2.04 for 
the three different times. The lengths of arrows are are nor- 
malized to 0.3c (left) and 0.1c (right). The thick solid lines 
denote the location of the apparent horizon. 
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FIG. 7. Fraction of the rest mass inside a coordinate radius 
r as a function of t/P for star (D) with AK/K = 1% (evolved 
with type I modified AMD gauge). R e denotes the initial 
coordinate length of the semi major axis. 



0.6 
0.5 
0.4 
0.3 



0.5 



1.5 



I I I 1 I I I I 1 I I I I 1 I , I I 1 I I I ,_ 






x— ••• ■ 


, , l\, , , 


1 1 1 1 1 1 1 


1 , , , 



2.5 




2.5 



FIG. 8. a and p at r = as a function of t/P for star 
(A) of various initial pressure depletion factors AK/K. The 
solid, dotted, dashed and dotted-dashed lines denote the cases 
where AK/K = 0, 1%, 5%, and 10%, respectively. 
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FIG. 9. Same as Fig. 6, but for star (A) with 
AK/K = 10%. The contour lines are drawn for 
p»/p, c = 1CT ' 3J for j = 0, 1, 2, • ■ • , 10 where p* c is 0.012, 
0.20, and 1.01 for the three different times. 
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FIG. 10. Same as Fig. 6, but for star (A) with 
AK/K = 10% and the nonaxisymmetric perturbation (3.2). 
The contour lines are drawn for c = W^°' 3j for 

j = 0, 1, 2, • • • , 10 where p* c is 0.012, 0.15, and 0.99 for the 
three different times. 
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FIG. 11. Same as Fig. 7, but for star (A), AK/K = 10%, 
with (solid line) and without (dashed line) the nonaxisym- 
metric perturbation (3.2). 



FIG. 13. Same as Fig. 7, but for collapse with the initial 
density profile (3.4) (mimicing star (A) driven into instability 
by accretion of additional matter). 
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FIG. 12. Same as Fig. 11, but for the mean square axial 
length (see Eq. (3.3). The solid and dashed lines denote sim- 
ulations with and without the nonaxisymmetric perturbation 
(3.2). 
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